PEC2 Análisis de datos ómicos
Máster interuniversitario de Bioestadística y Bioinformática
Información general sobre los datos
Durante el proceso de análisis, hemos encontrado que la base de datos correcta para realizar la anotación es la Ensembl en su versión 75 (2014), las versiones más recientes no contienen algunos de los genes que aparecen en el dataset. Por tanto utilizaremos esa para anotar los genes y calcular los GRanges para su ploteo posterior.
Selección aleatoria de los datos
En primer lugar vamos a seleccionar diez muestras aleatorias por cada grupo, primero seleccionaremos sobre el dataset de targets y luego cogeremos los conteos correspondientes a esas muestras.
# Only select and save the samples once, and set the sample seed just in case we have to resamples the same population
set.seed(28052020)
selected_data_file <- 'selected_data.Rdata'
if (!file.exists(selected_data_file)) {
targets <- read_csv('../data/targets.csv') %>%
group_by(Group) %>%
sample_n(10) %>%
arrange(Group) %>%
mutate(rn = Sample_Name) %>%
column_to_rownames('rn')
counts <- read_delim('../data/counts.csv', delim = ';') %>%
dplyr::select(X1, rownames(targets)) %>%
column_to_rownames('X1') %>%
as.matrix()
save(targets, counts, file = selected_data_file)
}else { load(selected_data_file) }
janitor::tabyl(targets, Group) %>% kable()| Group | n | percent |
|---|---|---|
| ELI | 10 | 0.3333333 |
| NIT | 10 | 0.3333333 |
| SFI | 10 | 0.3333333 |
Tabla 1: Muestras por grupo
Tabla 2: Muestras seleccionadasLlegados a este punto vemos que tenemos cargados el número correcto de muestras por grupo
Filtrado inicial
Antes de comenzar a explorar los datos vamos a eliminar los conteos muy bajos. Eliminaremos todos los genes que a lo largo de todas las muestras tengan menos de 10 observaciones. Aunque el sistema es un poco simple, la idea es simplemente reducir la carga computacional y eliminar el ruido más elemental, a la hora de analizar los datos, DSeq se encargará de eliminar los conteos más ruidosos y además utilizaremos log fold shrinkage para reportar los resultados y penalizar los datos más ruidosos.
Control de calidad y exploración
Vamos en primer lugar a controlar la calidad de los datos y a realizar una pequeña exploración de los mismos.
Antes de empezar vamos a cambiar la matriz de conteos y vamos a transformarla en formato largo, deberíamos de ser precavidos ya que esto funcionará bien en un dataset pequeño pero en dataset más grandes podría dar problemas.
counts_long <- counts %>%
as_tibble() %>%
rownames_to_column(var = 'gene_id') %>%
pivot_longer(-gene_id, names_to = 'Sample_Name', values_to = 'count') %>%
inner_join(dplyr::select(targets, Sample_Name, Group), by = 'Sample_Name') %>%
mutate(Group = factor(Group, levels = c('NIT', 'SFI', 'ELI'))) %>%
arrange(Group)Distribución de conteos por muestras y grupos.
Vamos a comprobar la distribución de conteos por muestra y grupo, primero estudiando su densidad y luego en forma de boxplot.
cowplot::plot_grid(
ggplot(counts_long, aes(x = log10(count + 1), color = Sample_Name, group = Sample_Name)) +
geom_density(alpha = .2) +
lmft +
theme(legend.position = 'None') +
labs(title = 'Conteo por muestra'),
ggplot(counts_long, aes(x = log10(count + 1), color = Sample_Name, group = Sample_Name)) +
geom_density(alpha = .2) +
facet_wrap(. ~ Group) +
lmft +
theme(legend.position = 'None') +
labs(title = 'Conteo por grupos'),
nrow = 2
)Figura 1: Plots de densidad de los valores de conteo por muestra y grupo
ggplot(counts_long, aes(y = log2(count + 1), x = glue::glue('{Group}_{Sample_Name}'), fill = Group, Group = Group)) +
geom_boxplot(alpha = .2) +
lmft +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(title = 'Conteo por muestra', x = 'Muestra') Figura 2: Boxplot de los valores de conteo por muestra y grupo
Vamos ahora a controlar el conteo total por cada muestra
ggplot(counts_long %>%
mutate(Sample_Name = glue::glue('{Group}_{Sample_Name}')) %>%
group_by(Sample_Name) %>%
summarise(count = sum(count)),
aes(x = Sample_Name, y = (count / 1e6))) +
geom_bar(stat = 'identity', alpha = 0, color = 'black') +
lmft +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) Figura 3: Histograma del número total de conteos por por muestra y grupo
Parece que tenemos una muestra del grupo ELI que está ligeramente desplazada respecto de las otras muestras y además tiene muchas menos medidas. En concreto es la muestra que termina en 2TC6J.
Como no tenemos acceso a los datos iniciales no podemos comprobar cual ha sido el fallo en la muestra, por tanto lo que vamos a hacer es eliminarla y obtener una muestra nueva de nuestra base de datos.
# Remove wrong sample
wrong_sample <- targets[str_detect(targets[['Sample_Name']], '2TC6J'), 'Sample_Name']
targets <- targets[!str_detect(targets[['Sample_Name']], '2TC6J'),]
# Replace sample
new_sample <- read_csv('../data/targets.csv') %>%
dplyr::filter(Group == 'ELI', Sample_Name != wrong_sample, !(Sample_Name %in% targets[['Sample_Name']])) %>%
sample_n(1)
targets <- bind_rows(targets, new_sample) %>%
mutate(rn = Sample_Name) %>%
column_to_rownames('rn')
# Recreate counts
counts <- read_delim('../data/counts.csv', delim = ';') %>%
dplyr::select(X1, rownames(targets)) %>%
column_to_rownames('X1') %>%
as.matrix()
# Refilter
counts <- filter_lc(counts, 1)Vamos a volver a hacer el control de calidad con histogramas
counts_long <- counts %>%
as_tibble() %>%
rownames_to_column(var = 'gene_id') %>%
pivot_longer(-gene_id, names_to = 'Sample_Name', values_to = 'count') %>%
inner_join(dplyr::select(targets, Sample_Name, Group), by = 'Sample_Name') %>%
mutate(Group = factor(Group, levels = c('NIT', 'SFI', 'ELI'))) %>%
arrange(Group)cowplot::plot_grid(
ggplot(counts_long, aes(x = log10(count + 1), color = Sample_Name, group = Sample_Name)) +
geom_density(alpha = .2) +
lmft +
theme(legend.position = 'None') +
labs(title = 'Conteo por muestra'),
ggplot(counts_long, aes(x = log10(count + 1), color = Sample_Name, group = Sample_Name)) +
geom_density(alpha = .2) +
facet_wrap(. ~ Group) +
lmft +
theme(legend.position = 'None') +
labs(title = 'Conteo por grupos'),
nrow = 1
) Figura 4: Plot de densidad de valor de conteos por muestra y grupo (Muestra Corregida)
ggplot(counts_long, aes(y = log2(count + 1), x = glue::glue('{Group}_{Sample_Name}'), fill = Group, Group = Group)) +
geom_boxplot(alpha = .2) +
lmft +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(title = 'Conteo por muestra', x = 'Muestra')Figura 5: Boxplot de conteos por muestra y grupo (Muestra Corregida)
ggplot(counts_long %>%
mutate(Sample_Name = glue::glue('{Group}_{Sample_Name}')) %>%
group_by(Sample_Name) %>%
summarise(count = sum(count)),
aes(x = Sample_Name, y = (count / 1e6))) +
geom_bar(stat = 'identity', alpha = 0, color = 'black') +
lmft +
theme(axis.text.x = element_text(angle = 45, hjust = 1))Figura 5: Boxplot de conteos por muestra y grupo (Muestra Corregida)
Parece que hemos corregido el problema que teníamos con nuestra muestra, por tanto proseguimos con el análisis. De ahora en adelante siempre utilizaremos la muestra corregida.
DSeqDataset
Llegado a este punto vamos a transformar el set de datos en un DSeqDataset. Antes de realizar este paso vamos a transformar el campo Group, el de interés para nuestro análisis, en un factor cuyo primer nivel será NIT que será el de referencia cuando analicemos los datos. Lo haremos en este punto para poder transformar los datos en el siguiente paso.
En este punto vamos a añadir los GRanges al objeto DSeq, es importante recordar que la base de datos correcta para realizar esta anotación es la v75, que es la que contiene los genes que aparecen analizados aquí.
targets[['Group']] <- factor(targets[['Group']], levels = c('NIT', 'SFI', 'ELI'))
dds <- DESeqDataSetFromMatrix(countData = counts, colData = targets, design = ~Group)
names_genes_dds <- str_replace_all(names(rowRanges(dds)), '\\..*', '')
genes_ensemb <- genes(EnsDb.Hsapiens.v75)
rowRanges(dds) <- genes_ensemb[names_genes_dds,]Transformación
Algunos de los métodos de exploración que vamos a utilizar funcionan mejor bajo condiciones de homocedasticidad. Sin embargo esto no es completamente cierto para nuestros datos ya que la varianza aumenta en función del conteo por gen, es decir los genes con mayor conteo tendrán mayor variabilidad.
Podemos verlo claramente en nuestros datos.
ggplot(counts_long %>%
group_by(gene_id) %>%
summarise(sd_gene = sd(count), rank_gene = mean(count)) %>%
mutate(rank_gene = row_number(rank_gene))
, aes(x = rank_gene, y = sd_gene)) +
geom_point() +
ylim(0, 2^10)Figura 6: Scatter plot entre los genes ordenados de acuerdo a su conteo y su desviación standard. Nótese que el gráfico está cortado en su parte superior.
Como vemos en la Figura 6la variabilidad aumenta con la media.
Vamos por tanto a transformar nuestros datos para reducir la heterocedasticidad. Como son pocas muestras utilizaremos variance stabilizing transformation o vst. Lo haremos con la opción blind que evitará que en el computo de la transformación influya la varianza debida a las condiciones experimentales.
Vamos a comprobar que hemos arreglado este problema y la variabilidad de los datos no depende, o lo hace menos de la media.
ggplot(assay(vsd) %>%
as_tibble(rownames = 'gene_id') %>%
pivot_longer(-gene_id, values_to = 'count', names_to = 'Sample_Name') %>%
group_by(gene_id) %>%
summarise(sd_gene = sd(count), rank_gene = mean(count)) %>%
mutate(rank_gene = row_number(rank_gene))
, aes(x = rank_gene, y = sd_gene)) +
geom_point()Figura 7: Scatter plot entre los genes ordenados de acuerdo a su conteo y su desviación standard. Nótese que el gráfico al contrario que el anterior, no está cortado en su parte superior.
Como vemos en la Figura 7, aunque no es completamente perfecto, la variabilidad esá menos asociada al conteo por gen. Vamos ahora a explorar los datos.
Distancia entre muestras
sampleDists <- dist(t(assay(vsd)))
library("pheatmap")
library("RColorBrewer")
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- vsd$Group
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette(rev(brewer.pal(9, "Blues")))(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows = sampleDists,
clustering_distance_cols = sampleDists,
col = colors)Figura 8: Plot de distancia entre genes un clustering jerárquico.
El agrupamiento es más fuerte entre las muestras del grupo ELI mientras que el resto no parecen agruparse particularmente bien.
PCA
Figura 9: Análisis de componentes principales. Los colores representan los diferentes grupos experimentales.
En el caso del PCA parece más prometedor, existe una diferencia clara entre condiciones en el eje X y además la cantidad de varianza explicada por ese factor es elevada, un 59%.
MDS
mds <- as.data.frame(colData(vsd)) %>%
cbind(cmdscale(sampleDistMatrix))
ggplot(mds, aes(x = `1`, y = `2`, color = Group)) +
geom_point(size = 3) +
coord_fixed()Figura 10: Plot de Multidimensional Scaling.
Como en el caso de las distancias parece que el grupo ELI se diferencia claramente de los otros mientras que SFI y NIT se mezclan más. Esto es coherente ya que Multidimensional scaling intenta situar los puntos en el plano de tal manera que tengan una distancia igual a la matriz de distancias que hemos calculado previamente.
Agrupamiento de genes
Vamos a seleccionar los 20 genes más variables a través de las muestras de los datos transformados sobre vst
topVarGenes <- head(order(rowVars(assay(vsd)), decreasing = TRUE), 20)
mat <- assay(vsd)[topVarGenes,]
mat <- mat - rowMeans(mat)
anno <- as.data.frame(colData(vsd)[, c("Group", 'Group')])
pheatmap(mat, annotation_col = anno)Figura 11: Plot de variabilidad entre los 20 genes con clustering jerárquico.
El único cluster que parece asociado a las condiciones está en la esquina superior izquierda. Parece que contiene un set de genes asociados a la condición ELI. EL cluster central no parece particularmente asociado a ninguna condición.
Conclusiones
Parece que los datos son correctos y la exploración de los datos también indica que podemos proceder con el análisis de expresión diferencial.
Análisis de expresión diferencial (DSeq) y Anotación
Como vamos a aplicar DSeq no es necesario aplicar ningún tipo de transformación por lo tanto trabajaremos directemente sobre los conteos, en el objeto DSeqDataset que hemos creado en el apartado anterior. En este caso no hemos añadido la variable sex al modelo, sin embargo sería correcto comprobar contra un modelo con y sin esta variable para saber si es necesario introducirla en el modelo, o quizás incluirla de todos modos ya que la influencia del sexo es probablemente cierta, se podría consultar con el investigador. Como la PEC no indica
Si quisiésemos por ejemplo trabajar con limma deberíamos de transformar los datos (VOOM) ya que no está preparada para conteos.
En este caso disponemos de tres condiciones:
- NTI: Non-Infiltrated Tissues
- SFI: Small Focal Infiltrates
- ELI: Extensive Lymphoid Infiltrates
Vamos a aplicar varios enfoques en este caso:
Compararemos las dos condiciones experimentales (SFI y ELI) contra la condición control (NTI) y luego compararemos las condiciones experimentales entre ellas.
En todo los casos utilizaremos como corte de significación p<.05 y en los casos de FDR un 5%. En otro caso se indicará explícitamente.
dds_pair_file <- 'dds_pair.RData'
if (!file.exists(dds_pair_file)) {
dds_pair <- DESeq(dds, parallel = T)
save(dds_pair, file = dds_pair_file) }else {
load(dds_pair_file)
}Extraemos los resultados para cada contraste de interés y luego los plotearemos uno por uno. También vamos a calcular los resultados con un shrinkage de el fold que nos ayudará a corregir la variabilidad en los casos en el que conteo sea muy bajo o muy variable y podremos ordenar mejor los resultados.
En el caso de los contrastes contra el nivel de referencia utilizaremos apeglm para el shrinkage y ashr, para el caso del ELI y SFI.
Para realizar las anotaciones de los genes utilizaremos una vez mas la Ensembl v75 en lugar de la org.Hs.eg.db que hemos utilizado en la PEC anterior.
annot <- function(rs) {
if(class(rs)=='GRanges'){
k <- rownames(mcols(rs))
}else{
k <- rownames(rs)
}
rs$symbol <- mapIds(EnsDb.Hsapiens.v75,
keys = str_replace_all(k, '\\..*', ''),
column = "SYMBOL",
keytype = "GENEID",
multiVals = "first")
rs$entrez <- mapIds(EnsDb.Hsapiens.v75,
keys = str_replace_all(k, '\\..*', ''),
column = "ENTREZID",
keytype = "GENEID",
multiVals = "first")
return(rs)
}
rs_dds_pair_file <- 'rs_dds_pair.RData'
if (!file.exists(rs_dds_pair_file)) {
k <-rs_SFI_NIT <- annot(results(dds_pair, name = 'Group_SFI_vs_NIT', alpha = .05))
rs_ELI_NIT <- annot(results(dds_pair, name = 'Group_ELI_vs_NIT', alpha = .05))
rs_SFI_ELI <- annot(results(dds_pair, contrast = c('Group', 'SFI', 'ELI'), alpha = .05))
rs_SFI_NIT_LFC <- annot(lfcShrink(dds_pair, coef = 'Group_SFI_vs_NIT', type = 'apeglm'))
rs_ELI_NIT_LFC <- annot(lfcShrink(dds_pair, coef = 'Group_ELI_vs_NIT', type = 'apeglm'))
rs_SFI_ELI_LFC <- annot(lfcShrink(dds_pair, contrast = c('Group', 'SFI', 'ELI'), alpha = .05, type = 'ashr'))
rs <- bind_rows(mutate(as_tibble(rs_SFI_NIT, rownames = 'gene'), contrast = 'Group_SFI_vs_NIT'),
mutate(as_tibble(rs_ELI_NIT, rownames = 'gene'), contrast = 'Group_ELI_vs_NIT'),
mutate(as_tibble(rs_SFI_ELI, rownames = 'gene'), contrast = 'Group_SFI_vs_ELI'))
rs_LFC <- bind_rows(mutate(as_tibble(rs_SFI_NIT_LFC, rownames = 'gene'), contrast = 'Group_SFI_vs_NIT'),
mutate(as_tibble(rs_ELI_NIT_LFC, rownames = 'gene'), contrast = 'Group_ELI_vs_NIT'),
mutate(as_tibble(rs_SFI_ELI_LFC, rownames = 'gene'), contrast = 'Group_SFI_vs_ELI'))
save(rs_SFI_NIT, rs_ELI_NIT, rs_SFI_ELI,
rs_SFI_NIT_LFC, rs_ELI_NIT_LFC, rs_SFI_ELI_LFC,
rs, rs_LFC,
file = rs_dds_pair_file) }else {
load(rs_dds_pair_file)
}Hemos calculado también los resultados con log fold shrinkage para observar como determinados genes con grandes cambios en el log fold van asociados a bajos conteos y alta variabilidad. LFCShrink penaliza estos genes reduciendo su log fold. Utilizaremos estos últimos para reportar los resultados. Además también hemos anotado todos los genes
Resultados
Annotations
Como hemos visto que hay muchos genes que no tienen anotación en org.Hs.eg, vamos primero a hacer un conteo rápido para saber cuantos genes mapeados tienen un correspondiente Symbol o EntrezID.
rs_SFI_NIT %>%
as_tibble(rownames = 'gene_id') %>%
transmute('Symbol' = ifelse(!is.na(rs_SFI_NIT[['symbol']])|rs_SFI_NIT[['symbol']]=='', 'Yes', 'No'),
'Symbol' = factor(Symbol, levels = c('Yes', 'No')),
'Entrez' = ifelse(!is.na(rs_SFI_NIT[['entrez']]), 'Yes', 'No'),
'Entrez' = factor(Entrez, levels = c('Yes', 'No'))
) %>%
janitor::tabyl(Symbol, Entrez) %>%
janitor::adorn_totals(where = c('row', 'col')) %>%
janitor::adorn_title() Entrez
Symbol Yes No Total
Yes 22762 20649 43411
No 0 0 0
Total 22762 20649 43411
Tabla 3: Tabla de éxito de mapeado de Ensembl a EntrezID y Symbol
Como vemos aproximadamente la mitad de los genes no tienen ‘EntrezID’ de acuerdo a la base de datos mientras que todos tienen un Symbol asociado.
Tablas de resultados
Vamos a presentar los resultados en tablas que son ordenables y filtrables para cada uno de los contrastes por separado. Hemos eliminado las filas cuyos p.valores correspondían a NA. Los datos de cada contraste se dividen en Log Fold Positivo y Negativo.
table_report <- function(df) {
one_table <- function(df, fun) {
df %>%
as_tibble(rownames = 'gene_id') %>%
fun() %>%
arrange(padj) %>%
drop_na(-symbol, -entrez) %>%
sample_n(100) %>%
mutate_if(is.numeric, function(x) { signif(x, digits = 3) })
}
return(list(one_table(df, fun = function(x) { dplyr::filter(x, log2FoldChange > 0) }),
one_table(df, fun = function(x) { dplyr::filter(x, log2FoldChange < 0) })))
}
short_sum <- function(df, alpha = .05) {
df %>%
as_tibble(rownames = 'gene_id') %>%
transmute(Regulated = if_else(log2FoldChange > 0, 'Up', 'Down'),
Regulated = factor(Regulated, levels = c('Up', 'Down')),
Significant = if_else(!is.na(padj),
if_else(padj < alpha, 'Yes', 'No'),
'Low Count'),
Significant = factor(Significant, levels = c('Yes', 'No', 'Low Count'))
) %>%
janitor::tabyl(Regulated, Significant) %>%
janitor::adorn_totals(where = c('row', 'col')) %>%
janitor::adorn_title()
}
print_datatable <- function(x, caption) { datatable(x %>% arrange(padj), filter = 'top', caption = caption) }Tabla SFI vs. NIT
| Significant | ||||
|---|---|---|---|---|
| Regulated | Yes | No | Low Count | Total |
| Up | 181 | 15330 | 7403 | 22914 |
| Down | 3 | 13585 | 6909 | 20497 |
| Total | 184 | 28915 | 14312 | 43411 |
Tabla ELI vs. NIT
| Significant | ||||
|---|---|---|---|---|
| Regulated | Yes | No | Low Count | Total |
| Up | 3403 | 13882 | 6643 | 23928 |
| Down | 1258 | 13081 | 5144 | 19483 |
| Total | 4661 | 26963 | 11787 | 43411 |
Tabla SFI vs. ELI
| Significant | ||||
|---|---|---|---|---|
| Regulated | Yes | No | Low Count | Total |
| Up | 1207 | 12892 | 6093 | 20192 |
| Down | 2817 | 13025 | 7377 | 23219 |
| Total | 4024 | 25917 | 13470 | 43411 |
Plots
Volcano plot
my_volcano <- function(df, alpha = .05, lim_fold = 1, title = '') {
df <- df %>%
tidylog::drop_na() %>%
mutate(sig_p = pvalue < alpha, sig_adj_p = padj < alpha) %>%
tidylog::drop_na()
unadj_p <- ggplot(df, aes(x = log2FoldChange, y = -log10(pvalue), color = sig_p)) +
geom_point() +
geom_hline(yintercept = -log10(alpha), linetype = 2) +
geom_vline(xintercept = lim_fold) +
geom_vline(xintercept = -lim_fold) +
labs(x = 'LogFold',
y = '-log10(p)',
title = paste(title, 'Sin corregir por FDR'),
color = paste0('p<', alpha)) +
facet_wrap(. ~ contrast) +
theme(legend.position = 'none') +
lmft
adj_p <- ggplot(df, aes(x = log2FoldChange, y = -log10(padj), color = sig_adj_p)) +
geom_point() +
ggthemes::theme_tufte(base_family = 'sans') +
geom_hline(yintercept = -log10(alpha), linetype = 2) +
geom_vline(xintercept = lim_fold) +
geom_vline(xintercept = -lim_fold) +
labs(x = 'LogFold',
y = '-log10(p_adj)',
color = paste0('p<', alpha),
title = paste(title, 'Corregido por FDR p-value')
) +
facet_wrap(. ~ contrast) +
lmft
return(list(unadj_plot = unadj_p, adj_plot = adj_p))
}cowplot::plot_grid(plotlist = my_volcano(df = rs, alpha = .05, lim_fold = 2, title = 'Sin LFC'), nrow = 2)cowplot::plot_grid(plotlist = my_volcano(df = rs_LFC, alpha = .05, lim_fold = 2, title = 'Con LFC'), nrow = 2)Figura 13: Volcano plot para cada una de las condiciones y los resultados con y sin shrinkage. Las líneas verticakes indican un logfold de 2 y la horizontal el límite de significación al 5%
La línea horizontal indica el límite de significación para un 5%, las líneas verticales un fold de 2, los colores indican si un gen ha resultado significativo o no.
Como vemos en el caso sin shrinkage aparecen multitud de genes con un fold muy elevado, sin embargo al aplicar LFC el número se reduce drásticamente, probablemente sea indicativo de abundantes genes con variabilidad elevada. La única comparación que no sale fuertemente penalizada es SFI vs. ELI. En la sección anterior hemos visto los volcano plots que hemos utilizado en primer lugar para comparar la diferencia entre los resultados cuando utilizamos un shrinkage y cuando no lo utilizamos. Por tanto no los repetiremos aquí
MA-Plot
Hemos ploteado para todos los contrastes tanto los resultados sin como con shrinkage. En el caso que no hemos un shrinkage podemos comprobar como la estimación del log fold es más variable con conteos bajos (logFold altos pero sin embargo no son significativos) y sin embargo se vuelve menos variable cuando aumenta el conteo. Una consecuencia de esto es que logFold de una determinada magnitud no son significativos si el conteo es bajo y sin embargo si lo son con conteos alto. En rojo aparecen coloreados los resultados significativos. Vemos también como la transformación vst es efectiva ya que reduce el logFold en los conteos bajos, con lo que en esa franja de los gráficos se concentran en cero.
SFI vs. NIT
Figura 14: MA plot de los resultados con y sin shrinkage
ELI vs. NIT
Figura 15: MA plot de los resultados con y sin shrinkage
Plot en el espacio genómico
Como ejemplo de gráfico en el espacio genómico vamos a mostrar el gen más significativo en la comparación SFI vs. ELI. Para ello primero reanalizamos los datos pero lo hacemos en formato GRanges y sacamos los fold changes de los resultados con shrinkage.
Vamos a plotear las 100,000 base alrededor del gen más significativo.
Utilizaremos tracks para mostrar: - Los p.valores ajustados (-log10) - Los log2fold corregidos - Los genes significativos al 5% - Los genes cuyo log2fold absoluto corregido es mayor que 2
resGR <- annot(results(dds_pair, name = "Group_SFI_vs_NIT", format = "GRanges"))
resGR$log2FoldChange <- rs_SFI_ELI_LFC$log2FoldChange
resGR$logpadj <- -log10(resGR$padj)
window <- trim(resGR[which.min(resGR$padj)] + 1e5)
strand(window) <- '*'
resGRsub <- resGR[resGR %over% window]
naOrDup <- is.na(resGRsub$symbol) | duplicated(resGRsub$symbol)
resGRsub$group <- ifelse(naOrDup, names(resGRsub), resGRsub$symbol)
status_p <- factor(ifelse(resGRsub$padj < 0.05 & !is.na(resGRsub$padj), "sig", "notsig"))
status_2 <- factor(ifelse(abs(resGRsub$log2FoldChange) > 2, "sig", "notsig"))
options(ucscChromosomeNames = FALSE)
g <- GenomeAxisTrack()
p_annot <- AnnotationTrack(resGRsub, name = "significant - gene ranges (padj<.05)", feature = status_p)
log_annot <- AnnotationTrack(resGRsub, name = "logfold > 2 - gene ranges", feature = status_2)
p_data <- DataTrack(resGRsub, data = "logpadj", baseline = 0,
type = "h", name = "-log10(padj)", strand = "+")
log_data <- DataTrack(resGRsub, data = "log2FoldChange", baseline = 0,
type = "p", name = "log2 fold change", strand = "+")
plotTracks(list(g, p_data, log_data, p_annot, log_annot), groupAnnotation = "group",
notsig = "deepskyblue", sig = "tan2")Figura 16: Plot track del gen más significativo y las 100,000 bases vecinas. En narajan aparecen los genes signficativos al 5% (track superior) o con un log fold superio a 2 (track inferior).
Podemos observar como en este área la mayoría de los genes tienen un log fold inferior a cero por tanto su conteo es menor en la condición ELI que en la condición SFI y existen varios genes significativos alrededor. También podemos observar como incluso genes con un fold change elevado no son significativos, probablemente debido a la variabilidad.
Análisis de enriquecimiento
Hemos utilizado el paquete clusterProfiler
- Para estos análisis hemos seleccionado todos aquellos genes cuyo p-valor corregido con FDR es menor a .1
- Hemos realizado un análisis por contraste
- Hemos utilizado las bases de datos KEGG y GO
- En el caso de KEGG hemos corregido siempre con un p<.05 y un q<.05 para ORA y p<.05 para GSEA
- En el caso de GO hemos corregido siempre con un p<.05 y un q<.05, y p<.05 para GSEA.
- Hemos realizado 10000 permutaciones en los GSE
Enlaces a los tipos de gráficos utilizados:
En este sección utilizaremos dos tipos de análisis sobre dos bases de datos diferentes:
- Overrepresentation y Gene Set Enrichment Analysis
Interpretación de los resultados
Para los resultados de cada uno de estos análisis presentamos el mismo conjunto de tablas y resultados. Una tabla que muestra los resultados obtenidos, con información relevante sobre cada uno de ellos y un enlace a cada uno de los elementos de la GO o KEGG encontrados. Los tipos de gráficos varían según el tipo de análisis.
En el caso de KEGG cuando una ruta resulta significativa, el diagrama de la ruta se descarga y se colorea dependiendo de la expresión de los genes incluidos en la ruta para el análisis.
enrich_set <- rs_LFC %>%
drop_na() %>%
dplyr::select(p.fdr = padj, ENTREZID = entrez, fold = log2FoldChange, contrast) %>%
dplyr::filter(p.fdr < .05)
to_entrezid <- function(df) {
df %>% dplyr::pull(ENTREZID)
}
to_fold <- function(df) {
sort(tibble::deframe(df %>% dplyr::select(ENTREZID, fold)), decreasing = T)
}
enrich_kegg <- function(cont, gene_list, pvaluecutoff = .05, qvaluecutoff = .05) {
gene_list_contrast <- to_entrezid(gene_list %>% dplyr::filter(contrast == cont))
gene_list_kegg <- bitr_kegg(gene_list_contrast,
fromType = 'ncbi-geneid',
toType = 'kegg', organism = 'hsa')$kegg
enrich_res <- enrichKEGG(gene_list_kegg,
organism = 'hsa',
pvalueCutoff = pvaluecutoff,
keyType = 'kegg',
qvalueCutoff = qvaluecutoff)
# Ok con el ncbi-geneid! https://www.genome.jp/dbget-bin/www_bget?dme:Dmel_CG14941 Ejemplo
return(enrich_res)
}
cont <- 'Group_SFI_vs_NIT'
gene_list <- enrich_set
gse_kegg <- function(cont, gene_list, pvaluecutoff = .05, n_perm = 10000) {
gene_list_contrast <- to_fold(gene_list %>% dplyr::filter(contrast == cont))
names(gene_list_contrast) <- bitr_kegg(names(gene_list_contrast),
fromType = 'ncbi-geneid',
toType = 'kegg', organism = 'hsa')$kegg
enrich_res <- tryCatch({ gseKEGG(gene_list_contrast,
organism = 'hsa',
pvalueCutoff = pvaluecutoff,
nPerm = n_perm,
keyType = 'kegg') },
error = function(cond) { tibble(ID = character(),
Description = character(),
SetSize = numeric(),
enrichmentScore = numeric(),
NES = numeric(),
pvalue = numeric(),
p.adjsut = numeric(),
qvalues = numeric(),
rank = numeric(),
leading_edge = character()
) })
# Ok con el ncbi-geneid! https://www.genome.jp/dbget-bin/www_bget?dme:Dmel_CG14941 Ejemplo
return(enrich_res)
}
enrich_go <- function(cont, gene_list, pvaluecutoff = .05, qvaluecutoff = .05) {
gene_list_contrast <- to_entrezid(gene_list %>% dplyr::filter(contrast == cont))
enrich_res <- enrichGO(gene = gene_list_contrast,
OrgDb = org.Hs.eg.db,
keyType = 'ENTREZID',
ont = "CC",
pAdjustMethod = "BH",
pvalueCutoff = pvaluecutoff,
qvalueCutoff = qvaluecutoff)
enrich_res <- setReadable(enrich_res, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
return(enrich_res)
}
cont <- 'Group_SFI_vs_NIT'
gse_go <- function(cont, gene_list, n_perm = 10000, pvaluecutoff = 0.05) {
gene_list_contrast <- to_fold(gene_list %>% dplyr::filter(contrast == cont))
enrich_res <- tryCatch({ gseGO(geneList = gene_list_contrast,
OrgDb = org.Hs.eg.db,
keyType = 'ENTREZID',
ont = "CC",
nPerm = n_perm,
minGSSize = 100,
maxGSSize = 500,
pvalueCutoff = pvaluecutoff,
verbose = FALSE) },
error = function(cond) { tibble(ID = character(),
Description = character(),
GeneRatio = character(),
BgRatio = character(),
pvalue = numeric(),
p.adjust = numeric(),
q.value = numeric(),
Count = numeric()) }
)
if (nrow(enrich_res) > 0) {
enrich_res <- setReadable(enrich_res, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
}
return(enrich_res)
}plot_enrich_ORA <- function(enrich_df, contrast, type) {
if (nrow(enrich_df) == 0) {
return(invisible(NA))
}
p1 <- heatplot(enrich_df) + labs(title = paste(contrast))
p2 <- dotplot(enrich_df)
plot(p1)
plot(p2)
if (type == 'GO' && nrow(enrich_go_SFI) > 1) {
p3 <- goplot(enrich_df)
plot(p3)
# cowplot::plot_grid(p1,p2,p3, labels = 'AUTO', nrow = 3)
}
}
plot_enrich_GSE <- function(enrich_df, contrast, type) {
if (nrow(enrich_df) == 0) {
return(invisible(NA))
}
if (type == 'GO') {
p1 <- emapplot(enrich_df) + labs(title = paste(contrast))
}else {
p1 <- ridgeplot(enrich_df) + labs(title = paste(contrast))
}
p2 <- dotplot(enrich_df) + labs(title = paste(contrast))
# cowplot::plot_grid(p1, p2, labels = 'AUTO', nrow=2)
plot(p1)
plot(p2)
}
enrich_go_table <- function(enrich_go) {
suppressWarnings(enrich_go %>%
tibble::as_tibble() %>%
dplyr::mutate(ID = glue::glue("<a href='http://amigo.geneontology.org/amigo/term/{ID}' target='_blank'>{ID}</a>", ID = ID)) %>%
dplyr::mutate_if(is.numeric, signif, digits = 2) %>%
dplyr::select(-one_of("geneID", "core_enrichment")) %>%
DT::datatable(escape = FALSE))
}
enrich_kegg_table <- function(enrich_kegg) {
if (nrow(enrich_kegg) != 0) {
enrich_kegg_tibble <- enrich_kegg %>% tibble::as_tibble()
gen_url <- function(x, pathID) { paste0("http://www.kegg.jp/kegg-bin/show_pathway?",
pathID, "/", x[pathID, "geneID"]) }
suppressWarnings(enrich_kegg_tibble %>%
dplyr::mutate(Link = lapply(ID, function(x) { stringr::str_sub(gen_url(x = enrich_kegg, pathID = x), 1, -2) }),
Description = paste0("<a href='", Link, "' target='_blank'>", Description, "</a>")) %>%
dplyr::mutate_if(is.numeric, signif, digits = 2) %>%
dplyr::select(-Link, -one_of("geneID", "core_enrichment")) %>%
DT::datatable(escape = FALSE)) }else
{ tibble() %>% DT::datatable(escape = FALSE) }
}
enrich_kegg_pathview <- function(enrich_kegg, gene_list, dir_path = './kegg_pathview') {
enrich_kegg_tibble <- enrich_kegg %>% tibble::as_tibble()
if (nrow(enrich_kegg_tibble) == 0) { return(invisible(NULL)) }
dev.null <- lapply(enrich_kegg_tibble$ID, function(x) { if (!file.exists(file.path(dir_path, paste0(x, '.xml')))) {
pathview(gene.data = to_fold(gene_list),
pathway.id = x,
species = 'hsa',
kegg.dir = dir_path
)
file.copy(from = paste0(x, ".pathview.png"),
to = file.path(dir_path, paste0(x, ".pathview.png")))
file.remove(paste0(x, ".pathview.png"))
} })
dev.null <- lapply(enrich_kegg_tibble$ID, function(x) { plot(cowplot::ggdraw() + cowplot::draw_image(image_trim(image_read(file.path(dir_path, paste0(x, ".pathview.png")))))) })
}Gene ontology
Over-representation analysis
enrich_go_SFI <- suppressMessages(suppressWarnings(enrich_go('Group_SFI_vs_NIT', enrich_set)))
enrich_go_ELI <- suppressMessages(suppressWarnings(enrich_go('Group_ELI_vs_NIT', enrich_set)))
enrich_go_SFI_ELI <- suppressMessages(suppressWarnings(enrich_go('Group_SFI_vs_ELI', enrich_set)))SFI vs. NIT
#### ELI vs. NIT
Gene Set enrichment analysis
gse_go_SFI_ELI <- suppressMessages(suppressWarnings(gse_go('Group_SIF_vs_ELI', enrich_set)))
gse_go_ELI <- suppressMessages(suppressWarnings(gse_go('Group_ELI_vs_NIT', enrich_set)))
gse_go_SFI <- suppressMessages(suppressWarnings(gse_go('Group_SFI_vs_NIT', enrich_set)))SFI vs. NIT
ELI vs. NIT
Resultados de GSE para GO
Resultados de GSE para GO
KEGG
En el caso de KEGG se han ocultado los pathways yo que en algunos casos son demasiados y hacen dificil navegar el documento. Se pueden desplegar en cada sección haciendo click en los puntos indicados.
Over-representation analysis
enrich_kegg_SFI <- suppressMessages(suppressWarnings(enrich_kegg('Group_SFI_vs_NIT', enrich_set)))
enrich_kegg_ELI <- suppressMessages(suppressWarnings(enrich_kegg('Group_ELI_vs_NIT', enrich_set)))
enrich_kegg_SFI_ELI <- suppressMessages(suppressWarnings(enrich_kegg('Group_SFI_vs_ELI', enrich_set)))SFI vs. NIT
ELI vs. NIT
Gene Set enrichment analysis
gse_kegg_SFI_ELI <- suppressMessages(suppressWarnings(gse_kegg('Group_SIF_vs_ELI', enrich_set)))
gse_kegg_ELI <- suppressMessages(suppressWarnings(gse_kegg('Group_ELI_vs_NIT', enrich_set)))
gse_kegg_SFI <- suppressMessages(suppressWarnings(gse_kegg('Group_SFI_vs_NIT', enrich_set)))SFI vs. NIT
ELI vs. NIT
Efectos ocultos
En este caso tenemos una variable oculta bastante evidente que es el sexo, vamos a comprobar si somos capaces de encontrarla con SVA
library("sva")
dat <- counts(dds_pair, normalized = T)
idx <- rowMeans(dat) > 1
dat <- dat[idx, ]
mod <- model.matrix(~ Group, colData(dds_pair))
mod0 <- model.matrix(~ 1, colData(dds_pair))
svseq <- svaseq(dat, mod, mod0, n.sv = 2)Number of significant surrogate variables is: 2
Iteration (out of 5 ):1 2 3 4 5
par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
stripchart(svseq$sv[, i] ~ dds_pair$sex, vertical = TRUE, main = paste0("SV", i))
abline(h = 0)
} Figura 17: Plot de las variables ocultas, en el eje x la variable real sex que no estaba incluida en el modelo.
No parece que la separación sea completamente perfecta, pero si parece que hay una separación aproximada particularmente en la primera variable oculta.
Archivos generados
Los archivos aparecen comprimidos en el repositorio en formato zip.
save_dose_result = function(enrich_df, fname){
enrich_df %>% tibble::as_tibble() %>% readr::write_csv(x = ., path = paste0(fname,'.csv'))
}
if(!file.exists('../results.zip')){
# Tables
rs_SFI_NIT %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_SFI_NIT.csv')
rs_ELI_NIT %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_ELI_NIT.csv')
rs_SFI_ELI %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_SFI_ELI.csv')
rs_SFI_NIT_LFC %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_SFI_NIT_LFC.csv')
rs_ELI_NIT_LFC %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_ELI_NIT_LFC.csv')
rs_SFI_ELI_LFC %>% as_tibble(rownames = 'gene_id') %>% write_csv('../results/results_SFI_ELI_LFC.csv')
# Enrichement
enrich_go_SFI %>% save_dose_result(fname='../results/ora_go_SFI')
enrich_go_ELI %>% save_dose_result(fname='../results/ora_go_ELI')
enrich_go_SFI_ELI %>% save_dose_result(fname='../results/ora_go_SFI_ELI')
gse_go_SFI %>% save_dose_result(fname='../results/gse_go_SFI_ELI_ELI')
gse_go_ELI %>% save_dose_result(fname='../results/gse_go_ELI')
gse_go_SFI_ELI %>% save_dose_result(fname='../results/gse_go_SFI_ELI')
enrich_kegg_SFI %>% save_dose_result(fname='../results/ora_kegg_SFI')
enrich_kegg_ELI %>% save_dose_result(fname='../results/ora_kegg_ELI')
enrich_kegg_SFI_ELI %>% save_dose_result(fname='../results/ora_kegg_SFI_ELI')
gse_kegg_SFI %>% save_dose_result(fname='../results/gse_kegg_SFI')
gse_kegg_ELI %>% save_dose_result(fname='../results/gse_kegg_ELI')
gse_kegg_SFI_ELI %>% save_dose_result(fname='../results/gse_kegg_SFI_ELI')
zip::zipr(files='../results',zipfile = '../results.zip')
}| Nombre |
|---|
| results_SFI_NIT.csv |
| results_ELI_NIT.csv |
| results_SFI_ELI.csv |
| results_SFI_NIT_LFC.csv |
| results_ELI_NIT_LFC.csv |
| results_SFI_ELI_LFC.csv |
| ora_go_SFI.csv |
| ora_go_ELI.csv |
| ora_go_SFI_ELI.csv |
| gse_go_SFI_ELI_ELI.csv |
| gse_go_ELI.csv |
| gse_go_SFI_ELI.csv |
| ora_kegg_SFI.csv |
| ora_kegg_ELI.csv |
| ora_kegg_SFI_ELI.csv |
| gse_kegg_SFI.csv |
| gse_kegg_ELI.csv |
| gse_kegg_SFI_ELI.csv |
Apéndice de código
Todo el código generado aparece en el markdown y se puede consultar o ocultar